The impact of diffusion on confined oscillated bubbly fluid 
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We consider the dynamics of monodisperse bubbly fluid confined by two plane solid walls and sub- 
jected to small-amplitude high-frequency transversal oscillations. The frequency these oscillations is 
assumed to be high in comparison with typical relaxation times for a single bubble, but comparable 
with the eigenfrequency of volume oscillations. A time-averaged description accounting for mutual 
coupling of the phases and the diffusivity of bubbles is applied. We find nonuniform steady states 
with the liquid quiescent on average. At relatively low frequencies accumulation of bubbles either 
at the walls or in planes oriented parallel to the walls is detected. These one-dimensional states are 
shown to be unstable. At relatively high frequencies the bubbles accumulate at the central plane 
and the solution is stable. 
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I. INTRODUCTION 

The dynamics of a single and multiple inclusions sus- 
pended in liquid medium has been attracting much atten- 
tion for many years. Of special interest is bubbly media 
with bubbles as generally soft, deformable objects. Be- 
cause of compressibility, bubbles are able to exhibit an 
additional "degree of freedom" if compared with solid, 
nondeformable inclusions. A simple example of a system 
where this factor becomes of crucial importance is bubbly 
fluid under high frequency oscillations. 

A well-known observation is the appearance of an aver- 
aged force on a single bubble in fluid under the action of 
acoustic fieldjii^i^ For instance, the time-averaged force 
exerted on the bubble of radius R in the standing wave 
of pressure p = Po{z) cos tot is given by 



pu;2 (02-1) 



(1) 



which in a particular case of po{z) — Pq cos kz results in 
F& = , ° sin(2fcz)e, (2) 



puj^ (f72 - 1) 
with a nondimensional parameter 



1 



puj^R^ 



2a 



(3) 



where fl presents the ratio of the eigenfrequency of vol- 
ume oscillations'*''^ to the frequency lu of external driving. 
Here, k = lu/cq is the wave number, cq is the speed of 
sound in the liquid free of bubbles, Pg is the mean pres- 
sure in the bubble, p is the fluid density, a is the surface 
tension, 7 is the adiabatic exponent, and = (0, 0, 1). 
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As it follows from expression the bubble tends to 
the antinodes of the pressure wave at low frequency lo 
{ft > 1) and to the nodes at high frequency w (SI < 1). 
This generic behavior is known as the primary Bjerknes 
effect and the force as in Eq. (P) is referred to as the 
Bjerknes force. 

The simplest approach to the averaged description of 
bubbly fluid is to treat the bubbles in a superimposed 
acoustic field individually, independent of each otheri^ 
Each bubble in the field experiences the Bjerknes force. 
However, such description lacks for possible collective (or 
feedback) effects and may fail even for very small con- 
centration of bubbles. The point is that a collection of 
bubbles infiuences the ambient so that eventually both 
phases can be firmly coupled, which is essential for the 
correct description. For instance, the presence of small 
amount of bubbles is known to qualitatively change prop- 
agation of acoustic wave in liquid.'' In a situation where 
the size of a bubble is small compared with the acoustic 
wavelength, the scattering on a single bubble is typically 
weak. However, an ensemble of bubbles is able to signif- 
icantly scatter the wave, because the bubbles coherently 
change their volume. In other words, in a liquid con- 
taining bubbles the speed of sound Cb can become much 
smaller than cq . If the acoustic wavelength is larger than 
the characteristic length L of the system, cq ^ ojL, the 
pure liquid behaves as incompressible. At the same time, 
in the bubbly medium it may happen that Cf, ^ luL and 
therefore scattering effects become importantj^i^ 

The impact of feedback effects on the averaged dy- 
namics of bubbly fluid has been addressed by Kobelev 
and Ostrovsky.''- They analyzed a coupled problem of 
the averaged drift of bubbles and scattering of acous- 
tic wave by the bubbles. Such factors as polydispersity 
of the bubbly fluid, dissipation of bubble oscillations and 
collisions of bubbles were taken into account. As a result, 
a generalized model of bubbly fluid has been obtained. 
Not only do the bubbles follow the prescribed averaged 
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force, their motion modifies the acoustic field and, hence, 
changes the averaged force. Particularly, propagation of 
a traveling acoustic wave in semi-infinite liquid for two 
situation has been analyzed. A bubbly layer is either 
of finite thickness or occupies the whole domain. In the 
former case, the so-called effect of self-transparency has 
been found. 

This study has been followed by a number of partic- 
ular analyses based on similar approximations with ac- 
count for cavitation and diffusion of gas from the bubbles 
into the liquid J^iiSiii^ It has been shown that a spatially 
uniform state and a one-dimensional soliton-like solution 
turn out to be unstable. As a result of self-organization, 
an asymmetric state emerges 

An essential point behind these studie o^°i^^i^^'^'^ is the 
assumption of the liquid quiescent on average. Although 
this approximation may be justified in the cited works, it 
becomes inappropriate in a number of situations, as e.g. 
in the present paper. Generally, one should go for averag- 
ing the momentum equation for the liquid phase without 
compromise. The first step in this direction has been 
recently performed in Ref. |14j. Both dissipation of the 
volume oscillation of the bubbles and bubble collisions 
are neglected. The frequency of vibrations is assumed 
to be so small that the liquid remains incompressible, 
the compressibility of the medium is caused solely by the 
bubbles. 

An expression of the averaged volume force has been 
obtained for monodisperse bubbly fluid. The theory is 
applied to study evolution of the initially homogeneous 
bubbly fluid in a thin layer confined by solid walls and 
subjected to transversal oscillations. Bubbles either ac- 
cumulate in planes parallel to the walls or settle at the 
boundaries. This accumulation process leads to infinite 
growth of the concentration, which makes the description 
invalid at a certain time. A more realistic picture cor- 
responds to saturation caused by dissipative processes, 
which have been ignored up to now. 

In the present paper we overcome this difficulty in a 
similar way we already applied for an incompressible sus- 
pension in the field of external force. We introduce dif- 
fusivity of bubbles, which naturally prevents the unphys- 
ical growth of the concentration and allows us to make 
a step beyond the previous findings. We start with the 
problem formulation in Sec. |TT1 Section IIIII focuses on 
the analysis of quasi-equilibrium states. The problem of 
stability is addressed in Sec. llVl and the results are sum- 
marized in Sec. El 



II. PROBLEM FORMULATION 

Consider monodisperse bubbly fluid flUing the space 
between two solid parallel planes separated by a distance 
2h. The system is subjected to transversal harmonic os- 
cillations of an amplitude a and a frequency u>. To apply 
the averaged description developed before;^ a number 
of requirements is to be satisfied. More precisely, we fo- 



cus on a dilute bubbly fluid with the equilibrium radius 
of the bubble R ^ h. Despite the smallness of volume 
fraction of bubbles, (j) 1, we describe the bubbles in 
terms of a finite field $ = <j)h'^/R'^, which is for simplicity 
referred to as the concentration. We consider small am- 
plitude and high-frequency oscillations in the sense that 
ah ^ and cuR'^ ^ i/, where 1/ is the kinematic viscos- 
ity of the fluid. More exact conditions that allow to ne- 
glect dissipative processes for a single oscillating bubble 
are discussed, e.g., in Ref. |7|. As it is mentioned before, 
we are interested in the situation where the frequency of 
external driving is comparable with the eigenfrequency of 
the breathing mode. We choose the Cartesian reference 
frame with the origin located in the central plane of the 
layer. Axes x and y are aligned in the central plane and 
axis z is normal to the solid boundaries. 

It has been shown beforCfi^ that under the above con- 
ditions a peaking regime occurs: the bubbles accumulate 
at certain planes, z = const, where their concentration 
grows abruptly to infinity within a finite time. As was an- 
nounced in SeclH this unphysical growth can be remedied 
by introducing diffusion. The generalization of the aver- 
aged model for diffusive bubbles is straightforward. The 
principal point is that the presence of diffusion does not 
influence the pulsation problem and enters the averaged 
equations only. As a result, diffusion appears naturally 
in the flux of bubbles [see Eq. (|4c)) ]. as one intuitively 
expects. 

By measuring the length, time, velocity, and pressure 
in the scales of h, h^D~^, Dh'^, and pvDh^^, where D 
is the bubble diffusivity, we arrive at the dimensionless 
boundary value problem: 



9u 

at 



u-Vu = -Vp + 3Qs*a$V?A^ 



(4a) 



5$ ,. . 



div u 

z 



0, j = Ud$-V$, (4b) 

0, Ud = u + gsVV'. (4c) 
±1: u = 0, e, •j = 0.(4d) 



Here u and are the velocities of the fluid and bubbles, 
respectively, p is the renormalized pressure, and j is the 
flux of bubbles. 

The amplitude ip of the velocity potential of pulsation 
flow, which enters Eqs. (Pa|) and ([3c]), is determined by a 
boundary value problem: 



^ - 0, 



f^2- 1 

z = ±l: • Vi/' = 1- 



(5a) 
(5b) 



For the sake of brevity, hereafter ip is called velocity po- 
tential. 

Boundary value problem (H])-® is governed by dimen- 
sionless parameters 



3 



and rj, given by ([3]). Here {(p) denotes the mean con- 
centration of bubbles. The first parameter, Qg, stands 
for the intensity of external driving. Parameter S is the 
Schmidt number, which is the ratio of the characteristic 
diffusion time to the viscous time scale. For most prac- 
tically relevant situations S is high. The third parame- 
ter, <I>Q, stands for feedback, it presents a measure of how 
strongly the fluid motion is influenced by the bubbles (for 
a similar situation, see Ref.lill). Technically, this param- 
eter serves as a scaling factor, it defines dimensionless 
concentration of the bubbles so that the space-averaged 
field $ is normalized by unity. As introduced in Sec HI be- 
low we distinguish two opposite cases of low (fJ > 1) and 
high [Vl < 1) frequencies. It is important to note that 
this distinction is purely conventional and has no contra- 
diction with the high-frequency approximation accepted 
for the averaged description.^i^ Generally, any value of ft 
satisfies this approximation. 

III. QUASI-EQUILIBRIUM STATE 

We now perform a one-dimensional analysis of a sta- 
tionary solution, in which all the fields are functions of 
the z-coordinate only. Although the averaged fluid veloc- 
ity is vanishing, the pulsation velocity is nontrivial. For 
this reason, this solution can be referred to as a quasi- 
equilibrium state or simply a quasi-equilibrium. We note 
that the vibration "Bjerknes" force exerted on the bub- 
bles does not vanish. However, in contrast to the previ- 
ous nondiffusive studyr^ this force is now compensated 
by the diffusive flux so that the total bubble flux jo turns 
to zero. 

As a result, Eqs. (SI)-® are reduced and for the quasi- 
equlibrium state we obtain 

$[, = gs'folV'o)', (6a) 

< = -^1^*0^0, (6b) 

z = ±1 : Vo = 1- (6c) 

Here primes denote derivatives with respect to z. A closer 
look at Eqs. ([6]) allows us to figure out symmetry prop- 
erties of the solution. Potential tpo and concentration $0 
have to be an odd and an even functions of z, respec- 
tively. Hence, the boundary value problem ^ can be 
treated in a half of the domain, say < z < 1, with a 
boundary condition 

z = 0: ^0 = (7) 

and the impermeability condition (|6c|) at 2 = 1. 
Next, Eq. (pa|) is easily integrated to yield 

$o = Cexp(Qs^o), (8) 

where the constant C is defined by the requirement of 
mass conservation: 

= / exp (Qs^o) dz. (9) 
Jo 



Note that from Eq. ([8]) and symmetry condition ([7]) it 
follows that the concentration has a maximum at 2 = 
for < 1 {Qs < 0) and a minimum in the opposite 
case, n > 1 {Qs > 0). This observation is in agreement 
with the well-known primary Bjerknes effect: bubbles 
accumulate in the nodes of pressure (or equivalently the 
velocity potential, as in our case) at high frequencies and 
in the antinodes at low frequencies. 

The substitution of solution ^ into Eq. ((6b)) leads to 
a nonlinear ordinary differential equation for the velocity 
potential: 

< = -§^,^'''^''^o. (10) 

Accounting for relation (|6ap. we integrate Eq. (fTO| and 
obtain 

where ipm = "00(1) and the result satisfies boundary con- 
dition ([5c)) . 

Equations (fTU)) and (fTT)) can be thought of as the sec- 
ond Newton law and the energy conservation law for a 
mechanical particle with ipo and z playing the role of the 
coordinate and time, respectively. This observation does 
not imply, however, full mechanical analogy because we 
deal with the boundary value problem but not the initial 
value problem as in mechanics. 

Although generally this problem can be solved only 
numerically, in a number of limiting cases we obtain an- 
alytical solutions. 

A. Low frequencies 

We first focus on the case of low frequencies, for which 
we introduce a parameter 



We start with the consideration of the limit of large 
Qs, in which all the bubbles accumulate at certain planes 
z = 2c or in other words form narrow bubbly screens. 
Outside these screens the fluid is almost free of bubbles. 
In such domains $0 = and therefore the Helmholtz 
equation (fTU)) [or Eq. (|6b)) ] is reduced to the Laplace equa- 
tion with a linear solution for ipQ. 

For < 2, the bubbles tend to the solid walls so 
that no bubbly screens appear away from the wall. The 
corresponding outer solution describing the potential in 
the bulk is 

V'^°)=/?z, f3 = -^ (12) 
1 — 

and consequently the concentration of bubbles is expo- 
nentially small. On the other hand, the inner solution, 
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which describes the bubbly screens locaHzed close to the 
walls, is given by the formulas: 



$n = 



?/3(l + /3) 

- msY 



MnF(C), 



F 



cosh/3^^ + /5"^sinh/3^^, 



(13) 

(14) 
(15) 



where ^ — Qs{l — z) is the "fast" coordinate near the 
wall. 

Thus, the concentration $o is high near the walls. In 
contrast to the case of nondiffusive bubblesji^ the bubbles 
now cannot leave the system. As a result, the dynamics 
of fluid is strongly influenced by the bubbles. The fluid 
moves as a solid body with the amplitude (3, which is 
larger than the amplitude of external driving. The "air 
cushions" formed of bubbles near the walls are akin to 
springs (see Fig. [1]) so that altogether the system acts 
as a resonator. Under periodic driving, the system dis- 
plays forced oscillations with the resonant value a — 1 
(/? oo), which separates two qualitatively different 
regimes. At values a < 1 {(3 > 0), the liquid at each 
point oscillates in phase with the walls. As it follows 
from (jlSp . function F(^) is monotonic for positive (3, and 
hence both the potential and bubble concentration are 
maximal directly at the walls. In the opposite case, a > 1 
(/? < 0), the liquid core oscillates in counter-phase with 
respect to the walls. Note that function F{£^) is no longer 
monotonic. Thus, although the concentration maximum 
is located very close to the wall, but not exactly at the 
wall. At the critical value a — 1, resonant amplification 
of oscillation occurs. In this particular situation, even 
small dissipation must be taken into account. 

For 2 < < 12, the bubbly screen is localized at 
the point z = zi = 2a^^. While far from this point the 
potential is a linear function 



\Z - Zl\ - Zl, 



(16) 



the solution close to the screen can be described as 

Qs 



$0 



cosh'^ S, ' 



ipQ = -zi + {ziQs) Mncosh^ 



(17a) 
(17b) 



with ^ = (z - Zl) ziQs- 



liquid core 



bubbles 



FIG. 1: Bubbly fluid as a resonator. At low frequencies bub- 
bles localize near the walls and become equivalent to springs, 
while the liquid plays the role of solid body. 



At larger the number of bubbly screens increases. 
For 2n(2n - 1) < < (2n -|- l)(2n -|- 2) there exist n 
bubbly screens localized at 

z = Zl = 2na^^, Z2 = 3zi, z„ — {2n — l)zi. (18) 

We note that in this case formulas p7|) remain valid 
in the vicinity of bubbly screen k {k — 1, . . . , n), with 
S^k — [z ~ Zk)ziQs and zi defined by Eq. (fTB]) . Another 
difference is that the sign of V'o*'' changes for even k. 

We now proceed to the opposite limiting case of small 
Qs, which is described by an asymptotic solution 



- ,/,(o) 



sinaz 



(1) 



*n = 1 



Ci = 



a cos a 
sin 2a — 2a 

4a'^ cos^ a 



(19a) 

Ci, (19b) 
{C^l + QsCi), (19c) 



where 



z cos az — (cos a — a sin a) '0g 



(0) 



1 



32q;'' cos^ a 



sin 3q!Z — 3a cos "iaij^^ 



(0) 



(20) 



with /o = {Aa^Ci cos^ a + 3)/(8a2 cos^ a). 

These results can be easily explained as follows. Small 
values of are equivalent to intensive diffusion. As 
a result, spatial inhomogeneities in the distribution of 
bubbles are smoothed out by diffusion. The Bjerknes 
force can lead to a small correction only, which results 
in a weakly nonuniform concentration field. Note that 
this quasi-equilibrium resembles the solution obtained for 
early stages of evolution in the nondiffusive approxima- 
tion [see formulas (71) and (73) in Ref. This similar- 
ity is caused by the initial conditions chosen in the form 
of uniformly distributed bubblesJ^ 

We next discuss numerical results. In Fig. [2] we show 
distributions of $o and -0o for a^ = 0.5. The dependen- 
cies are presented for different values of Qs- We note 
that the potential is linear everywhere, except for the 
vicinity of the wall. Because of the exponential depen- 
dence of $0 on ipQ, even a small change in profile ipoi^) 
drastically influences the concentration profile. To vali- 
date asymptotic solution for large Qs, we provide 
Fig. ^c). Here we demonstrate the variation of an aux- 
iliary field (f) = Qs^$o as a function of ^. It can be seen 
that even at Qs = 5 the numerical results are in good 
agreement with the asymptotic solution. 

Similar solutions are shown in Fig. [3] for = 1.7, 
when the maximum of concentration is located close to 
the wall, but not directly at it. This has been rigor- 
ously proved for large Qs- However, as it can be seen in 
Fig. [21Ia), a very similar situation takes place for finite 
values of Qs- For the values of Qs used in Fig. [3l the 
asymptotic solution is not as good as in the case in Fig. [21 
It should be emphasized, that the reliable agreement with 
the asymptotic solution is achieved at Qs > 50. 
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z z i 

FIG. 2: Quasi-equilibrium states at — 0.5. Profiles of the concentration of bubbles $o (a) and velocity potential t/jq (b) at 
Qs = 2, 5, 10, shown by dotted, dashed, and solid lines, respectively. Variation ol <}) = Qg^^o with ^ — Qs(l — z) for the same 
values of Qs and (c); circles present the asymptotic law according to formula ([T^. 
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FIG. 3: Quasi-equilibrium states at — 1.7. Profiles of the concentration of bubbles <l>o (a) and velocity potential ipo (b) at 
Qs = 2, 5, 10, shown by dotted, dashed, and solid lines, respectively. Variation of <^ = Qj^^o with ^ — Qs(l — z) for the same 
values of Qs and (c); circles present the asymptotic law according to formula (|13[) . 




FIG. 4: Maximal value of the bubble concentration, <I>m = 
maxz $o(z), as a function of Qg^. Solid lines show the nu- 
merical results, dashed lines are plotted according to formula 
(fT3l) . for Qs — > cxD. Lines 1 correspond to = 0.5, lines 2 - 
to = 1.7. 



The dependence of the concentration maximum on pa- 
rameter Qs is demonstrated in Fig.lH As before, one sees 
that the smaller is the value of a, the better asymptotic 
formula works. 

In Fig. [5] we plot numerically obtained profiles for 
larger values of . For — 4, 16, and 40 one, two, 
and three bubbly screens, respectively, exist in a half of 
the layer, < 2 < 1. The velocity potential is nearly a 
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FIG. 5: Profiles of the bubble concentration (a) and the po- 
tential of the pulsation velocity (b) plotted for Qs = 100. 
Solid, dashed, and dotted lines correspond to o? = 4, 16, 40, 
respectively. 



piecewise-linear function of z. 

B. High frequencies 

At high frequencies, 17 < 1, we introduce another aux- 
iliary parameter 



3$a 
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FIG. 6; Profiles of the bubble concentration and the velocity 
potential for different Qs and 6? . (a): The results correspond 
to 5? = 0.1, parameter Qs = —80 (dotted lines) and Qs ~ 
—280 (solid lines). On the scale of the figure, the solid and 
dotted lines for the potential xpo cannot be distinguished, (b): 
Similar dependencies for Qs = — 1. Parameter 5? = 10 (solid 
lines) and 5? = 20 (dotted lines). Asymptotical solution (|26|) 
valid for large is shown by circles. 



and recall that for high frequencies Qs < 0. 
In the limiting case \Qs\ = 3> 1 we obtain: 



^0 ^ ^0 ' ^0 



(1) - ;^,2fl,(o) 







/27r 



|exp(-e2^ (21a) 
(21b) 



il^Q = z + e^g{^), g= — r-erf^. 



(21c) 



where ^ = z/e and erf z = (2j^) exp{—y^)dy is the 
error function. This solution indicates that bubbles ac- 
cumulate at the center of the layer, z = 0, which corre- 
sponds to the node of the pulsation pressure. The veloc- 
ity potential is the same as for the pure liquid up to a 
small correction. For instance, the numerically obtained 
results shown in Fig. [DJa) perfectly match asymptotical 
solution (PTjl . 

In the opposite case, \Qs\ Ij the nonuniformity of 
concentration is small, $o is weakly enhanced at the cen- 
ter with a relative decrease near the boundaries. This 
case is described by an asymptotic solution of the form 



^0 



(0) 



',(1) 



sinh az 



8a^ cosh" a 
$0 = l + Qs (V'o + Ci) , Ci 



3 rV'o > Vo 
■y. 

la — sinh 2d 



a 



cosh( 



Ad^ cosh a 



,(22) 
(23) 



with 



(1) _ 



/o z cosh dz — (cosh d — d sinh d) tpQ 



(0) 



1 

4a 



sinh Zdz — 3a cosh idip^Q^ 



(24) 



where /o = (4a^Ci cosh^ a — 3). 

We point out an interesting case where the frequency of 
external driving only slightly exceeds the eigenfrequency 



of a single bubble, which corresponds to high values of 
a. Physically, this means that a boundary layer emerges 
near the wall. For this reason, we introduce the fast 
coordinate rj = a(l — z) and present the velocity potential 
as 



(25) 



Since potential is small, we notice from Eq. ([8]) that 
the concentration is nearly unity. This observation al- 
lows us to linearize Eq. (jlO|) and to figure out that 
/ = exp {—rj). As a result we obtain 



$0 



a 



§^ (6-3-7 „ 36-") 

8a3 



1 



Qs 



-2'n 



Qs 
253' 



(26a) 
(26b) 



A comparison of result (j26p with the numerical solu- 
tion is provided in Fig. [S^b). It is clearly seen that 
the asymptotical solution works well even at a = 20. 
With the increase of \Qs\ the difference between the an- 
alytical and numerical results becomes more pronounced. 
This tendency is easy to explain by looking at expression 
(|26bp . Larger values of \Qs\ mean stronger influence of 
the nonuniform part of concentration. 



IV. STABILITY ANALYSIS 

We now pose the question whether the solutions found 
in Sec. IIIII are stable. To answer this question, we intro- 
duce small perturbations of the bubble concentration 0, 
velocity potential of the pulsation motion \E', fluid veloc- 
ity U, and pressure P. By substituting the perturbed 
fields into Eqs. ^ and ([5]) and linearizing the problem 
with respect to small perturbations, we arrive at 



1 9U 
S~dt 

'dt 
F 



-vp-f-3$ar, divU==o, 



(27a) 



-U • V$o - div J, J = F - V0, (27b) 

Qs [2$oV (V^o*) + 0VV^'] , (27c) 
3$ 

" ($o* + ^o0), (27d) 

J = ■ V* = 0. (27e) 



fi2 _ 1 
±1 : U = e 



Since the base quasi-equilibrium state possesses O2 
symmetry, we do not have to treat the full three- 
dimensional problem. For this reason, we restrict our 
analysis by the two-dimensional stability problem. We 
assume that all the perturbation fields are independent 
of y and the corresponding component of the velocity 
vanishes, Uy = 0. As a result, for the two-dimensional 
incompressible velocity field we can introduce a stream- 
function if defined by relation 



U = V X (i^ej^ 



(28) 



where By = (0, 1, 0). 
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We apply operation Vx to Eq. (|27aP and consider the 
perturbations proportional to exjp[ikx + Xt). Here k is 
the real wavenumber and A is the complex growth rate. 
As a result, we obtain a boundary value problem for the 
z-dependent amplitudes of perturbations 

- i^V - 6iA:Q5$a*o ($0* - ^00) , (29a) 



S 



X(j) ^ -ifc$;,(^ - J' + (2Qs*o^o* - 0) (29b) 
DH = (^0^ + ^00), (29c) 



z = ±1: ip = ip' = J = = 



(29d) 



where = /dz^ — k"^ is the Fourier image of the 
Laplace operator and 

J EE 2Qs [$o (V'o^')' + V'oV'o'/'] - </>'• 

Having solved this boundary value problem, one finds the 
spectrum of eigenvalues A as a function of dimensionless 
parameters. The analytical solution can be obtained only 
in a few limiting cases, to solve the problem numerically 
we apply the standard shooting method. We note that 
in all our calculations A is found to be real. 

We emphasize that to consider a practically relevant 
limit of large Schmidt numbers, S ^ 1, the left hand 
side of Eq. (|29aP should be suppressed, which simplifies 
the analysis. Physically, this approximation implies a 
very fast relaxation of the perturbations associated with 
the flow. 



A. Low frequencies 

It can be easily shown that the quasi-equilibrium state 
is unstable for > 1 at arbitrarily small Qs- To prove 
this statement, let us have a look at the stability problem 
in the limit of small external driving, Qs «C 1, when the 
base state is defined by Eqs. p9)) . For the sake of brevity, 
below we omit superscript "(0)" for the leading part of 

the potential V'o"''- confusion is unlikely, because 

the first correction ip^^^ does not influence the further 
analysis. We expand the perturbations and the growth 
rate A in series with respect to 

(j) ^ (j)o + Qs0i + • • • , * = *o + Qs'^i + ■■■ (30a) 
(/? = V3o + Qsipi + ■ ■ . , (30b) 

and arrive to the zero order at a problem 

L^(^o = - A5-1) y'o = 0, (31a) 



)2 

,2 „2 



(31b) 

^'^ = {k'-a')^o, (31c) 
z = ±1 : ^0 = ¥^0 = (^0 = *o = 0. (31d) 

As we see, all the fields are decoupled, the solutions 
for ipo and 0o are given by sets of even and odd functions 



with real negative eigenvalues A. Note that the eigenval- 
ues associated with ifg are proportional to the Schmidt 
number. As for real bubbly fluids S is large, the pertur- 
bations of the flow decay extremely fast. The eigenvalue 
spectrum of the concentration, which is used in the fur- 
ther argumentation, are given by values 



A. 



0,1,2, 



(32) 



In other words, all the modes mentioned are decaying 
in time and therefore cannot lead to instability. Because 
we are now interested in growing and neutrally stable 
modes (A > 0), to this order we should put 



0. 



(33) 



The solution for the velocity potential is trivial, ^Pq = 
0, unless k — \J — ir'^m? /A, m = 0, 1, ... . We next 
deal with the simplest case of to = 0, which is the only 
option allowed for all possible a. In this case the bound- 
ary value problem for has a constant solution, which 
without loss of generality can be set to unity: 



*o = 1- 



(34) 



As we see, to the zero order no instability is detected 
and we proceed to the next order. Thus, in addition to 
expansions pop we present the wavenumber as 



k = a + Qski + . 
and obtain to the first order in Qs'. 



0, 



2aki ~ a $ 



2^(1) 



(35) 

(36a) 
(36b) 

(36c) 



z = ±1 : = = ^\ = 0, 0; = 2^-0. (36d) 

The solution for the streamfunction </ji = as before, 
whereas for the concentration of bubbles we obtain either 



2a2 + A 



2a2 + A 



2o?i)Q + A 



sinh qz 
q cosh q 



2a^tPo + A- 



sm qz 



q cos q 



A > -a^ 



A < 



(37) 



(38) 



where = + A and q^ = —o? — A. As it can be seen 
from relation (|32p for k = a, solution (|38p diverges at 
A ~ Xn for n odd. 

The solvability condition for Eq. (|36cp with appropri- 
ate boundary conditions can be obtained by integrating 
this equation across the layer, which for A > — yields 

, , o sin a cosh g — a cos a sinh (7 

fci = kg, + X (39) 

(2a^ + X) q cos a cosh q 
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FIG. 7: Stability diagram (a) for = 0.1 (lines 1) and = 
1 (lines 2). Solid and dashed lines correspond to numerical 
calculations and analytical result (|42p . respectively. Growth 
rates A as a function of fci (b). Dashed, solid, and dotted 
lines present the dependence at — 0.25, 1, 4, respectively. 
Inset provides a comparison of numerical results (circles) and 
analytical solution (|39} and (|40p (solid lines) for Qs = 0.2 
and — 1. 



and in the opposite case, A < —a^ we have 



ki — fcge ~ 



q sm a cos q — a cos a sm g 
(2(22 + A)^ 9 cos a cos g 



where 



2a — sin 2a 
^ 2 (2a2 + A) cos2 a 



> 0. 



(40) 



(41) 



The condition of neutral stability is defined by the re- 
quirement A = 0, which leads to ki = fcgeU^o = ^qe(O). 
Taking this observation into account in relation (j35p we 
figure out the border of stability to be 



(42) 



which is valid at small Qs and k ^ a. This result 
is in good agreement with numerical calculations, see 
Fig. [71 We indicate that the stability border is inde- 
pendent of both and S, even for finite Qs- As we 

(c) 

see, perturbations grow at any Qs > Qs ■ This growth 
takes place even for infinitely small intensity of external 
driving, where the perturbations are characterized by k 
slightly exceeding a. As a result, we conclude that at low 
frequencies, O > 1, the quasi-equilibrium state is always 
unstable. 

Let us come back to small values of Qs- We note that 
by setting A = in Eq. (|37p . one ends up with 0i = 2'ipQ. 
For the full concentration field we have [see relations (HH)] 



1 + Qs (^o' 



:Cexp \^Qs {^Pq + ^qY 
Ci) + 2Qs^o*o. 



(43) 



This fact indicates that for small Qs solution ([8]) remains 
valid even for the perturbed fields taken at the stability 
border, fci = fcge(O), with $ = ^o+Qsc/)! and V' = V'o+*o 
instead of $0 and ipoi respectively. Hence, the branching 



solution is another quasi-equilibrium state, but in con- 
trast to that in Sec. IIIII this state is two dimensional. 
This is a direct consequence of the specific form of func- 
tion $. For the concentration being an arbitrary func- 
tion of potential t/j, but the potential only, $ — F{ip), 
the feedback term in Eq. (|4ap can always be presented 
as gradient. Thus, this term redistributes pressure but 
does not generate the averaged fluid flow. Note that this 
result is valid even at finite values of fc — a and explains 
why the stability border is independent of S and $a for 
a fixed. These parameters enter Eqs. ^ and ([5]) only 
together with u. 

Let us now discuss the behavior of the growth rates as 
functions of /c. At /c w a these dependencies are described 
by Eqs. and (gOl), which are tabulated in Fig. El^b). 
The inset of this figure provides a comparison with the 
numerically obtained results. It can be seen that fci tends 
to infinity as A — > A„ for n odd, see relation ([5^ . This 
result is reasonable as it provides the matching of the 
different solutions separated by the critical value k — 
a- Next, it is clear from Fig. [Tfb) that in the vicinity 
of /c = a a rearrangement of branches occurs. Starting 
from A„ with n odd at fc > a, the growth rate steadily 
increases with the decrease of k and at A; < a reaches the 
value Xn-2- Moreover, a similar variation of the lowest 
odd branch, namely Ai, results in A +00 as A; —> a -1-0. 
Thus, the growth rate has a pole at fc = a and no positive 
growth rates exist in the spectrum at fc < a. For this 
reason, the domain with fc < a is marked as "stable" in 
Fig. Ha). 

This unstable mode originates from the problem of nat- 
ural oscillations for the velocity potential tfj (for the uni- 
form distribution of bubbles), $0 = 1- As we see from 
Eq. (|36b[) . this eigenmode induces the perturbations of 
concentration. Because of feedback, the concentration 
influences the potential and the system eventually be- 
comes unstable. This instability takes place for a base 
state with any nonvanishing ipo- The simplest example 
of such mode, inherent in Eqs. Q and is analyzed 
in Appendix [XI 



B. High frequencies 

We now consider the stability at high frequencies, fl < 
1. As before, in a few limiting cases, boundary value 
problem ([29[) admits an analytical solution. 

First, we focus on the limit of large \Qs\, when bubbles 
accumulate at the center of the layer and the potential 
of pulsation motion is nearly linear. This base state is 
described by Eqs. (|2ip . An accurate analysis of this situ- 
ation is performed by means of the matched expansions 
method (see Appendix [B[). which results in the spectrum 
of growth rates 



A„ = 2nQs, 71 = 0, 1, 2, . 



(44) 



Numerical results and asymptotic law agree well. 
The agreements becomes better for bigger n, see 
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FIG. 8; Growth rates at high frequencies plotted for 5? = 0.1, 
S = 100, $a = 1. (a): four lower branches of the spectrum at 
k = \, where solid lines present numerical results and dashed 
lines show the approximation for \Qs\ 3> 1, see formula (|44[) 
for n = 1, 2, 3. (b): Variation of the growth rate with k for 
Qs = —100 (solid line) and Qs ~ —280 (dotted lines); dashed 
line shows the asymptotic law according to (|45|l. 



FIG. 9: Growth rates at high frequencies presented for S = 
100, $a ~ 1, and k = 1. (a): Comparison of numerical data 
(solid line) and approximate formula (146 |l (dashed line), Qs = 
1. (b): Variation of the growth rates with Qs; parameter 
Se^ — 0.1, 1, 10 correspond to solid, dashed, and dotted lines, 
respectively. 



Fig. Ha). 

Except for n = 0, these branches display strong tem- 
poral decay of perturbations, which increases with the 
growth of driving intensity, \Qs\- Quite similar behavior 
of the spectrum has been recently observed for dielectric 
particles accumulated at the center of the layer under the 
action of dielectroporetic force^^^ 

For n = a more delicate analysis is needed. Refer- 
ring to Appendix [Bl for the details, we provide here the 
eventual result valid at the limit S ^ 1: 

, , , sinh^ k-k^ ^( 1 \ 

Ao = -k^ - 3$afc — 7—, r + O , . (45 

sinh2fc-2fc 

Figure O^b) shows the comparison of numerical results 
with approximation (|45|1 . Again, the results agree well, 
though with a slight distinction for higher k, where a 
correction to Aq becomes non- negligible. 

In another limiting case, ^ 1, when the base state 
is given by Eqs. (|26p . the largest growth rate is 

A.-.^(l + 2|f). ,46, 

Note that the Bjerknes force provides a small negative 
correction to the decay rate caused by difFusivity, so that 
the role of vibration force is destabilizing. As it becomes 
evident from Fig. [H^a) , formula works well even at 
a2 ^ 10. 

We have also checked several other limiting cases: 
$a ^ 1, when the there is no generation of the averaged 
flow, and ^ 1, when the potential of the pulsation mo- 
tion is linear. These analyses as well as numerical tests 
show that quasi-equilibrium state is stable. An example 
of calculations in which the Bjerknes force may become 
destabilizing is presented in Fig. [Hb). This destabiliza- 
tion, however, does not eventually lead to instability. 

Thus, our numerical and analytical results show that 
at high frequency the quasi-equilibrium state is stable. 



V. CONCLUSIONS 

We have considered the dynamics of monodisperse 
bubbly fluid confined by the plane solid walls. The 
system is subjected to small-amplitude high-frequency 
transversal oscillations. This frequency of external driv- 
ing is assumed to be high in comparison with typical 
relaxation times for a single bubble. At the same time, 
the ratio of the eigenfrequency of volume oscillations 
to the frequency of external driving, is of order unity. 
The time-averaged description developed in Ref. [H has 
been generalized. In contrast to the original model, we 
have taken into account the dijfusivity of bubbles, which 
allows us to prevent unbounded accumulation of bubbles 
found out earlier ."ii 

The quasi-equilibrium states, in which the fluid is qui- 
escent on average and the concentration of bubbles is 
nonuniform, have been systematically explored. In the 
state of quasi-equilibrium, the Bjerknes force, which acts 
on compressible bubbles, is balanced by the diffusive flux 
of bubbles. We stress that in contrast to the case of a 
single bubble, the ensemble of bubbles significantly in- 
fluences the characteristics of the fluid phase, which is 
referred to as feedback effects. Technically, this collective 
bubbly ensemble-induced effect is taken into account by 
coupling the phases without compromise. As a result, 
we are able to observe that the bubbles influence the 
pulsation field and therefore the Bjerknes force itself is 
changed. 

At a low frequency, Q > 1, we detect accumulation 
of bubbles either at the solid boundaries or in planes 
oriented parallel to the walls. Bubbly screens predicted 
in nondiffusive considerationjii are smeared by diffusion. 
As a result, the corresponding structures are stationary 
and no longer singular objects. We have shown that 
all these one-dimensional states turn out to be unsta- 
ble. What is interesting, the branching solution satisfies 
the criterium of quasi-equilibrium. This fact indicates 
that although the one-dimensional solutions are unsta- 
ble, two-dimensional quasi-equilibrium states and their 
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stability may become of interest. 

At a high frequency, fl < 1, the maximal value of the 
concentration is at the center plane of the system. As in 
the case of low frequencies, this peak can be very sharp, 
when the Bjerknes force dominates over the diffusive flux, 
or smooth in the opposite case. This one-dimensional 
state has been shown to be stable for any values of gov- 
erning parameters. 
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APPENDIX A: STABILITY OF UNIFORM 
OSCILLATIONS OF BUBBLY FLUID 

Consider motionless bubbly fluid, Uq = 0, which fills 
infinite space. We assume that bubbles arc uniformly 
distributed, $o = 1, and admit that ipo ^ 1. Recall 
that while obtaining the averaged model, the pressure 
pulsations were assumed proportional to the velocity po- 
tential, tp. Hence, physically, the assumption of t/jq — 1 
implies spatially uniform oscillations of the pressure field. 
We indicate that although such assumption is rather 
hypothetical from the practical point of view, it helps 
us to figure out the reason of the instability found in 
Sec. IIV Al Thus, in the system under consideration, the 
pressure oscillates with an amplitude H and frequency 
Lu, low in the sense f2 > 1. For this system, the pa- 
rameter characterizing the intensity of external driving 
is Qs = n2(2pw)-2[z.D(r!2 _ 

In order to investigate the stability of this state, we in- 
troduce small perturbations of the concentration, <j), and 
the potential of the pulsations, ^. After the linearization 
of Eqs. (HI) and (O with respect to the perturbations, one 
arrives at a problem 



ot 



0. 



(Ala) 
(Alb) 



We note that the perturbations of the flow effectively de- 
couple and turn out to decay. This is because for the 
case of interest the averaged vibration force in Eq. (Ha|) 
becomes gradient. Hence no averaged flow can be in- 
duced within the linear approximation. 



We seek the solution of Eqs. (jAip proportional to 
exp (Xt -|- ik • r) and obtain a dispersion relation 



A = -fc^ - 



2 2QsPa^ 



(A2) 



where k is the wavenumber. 

This relation qualitatively reproduces the picture of 
the instability shown in Fig. Eljb). As we can see, A is 
positive in a range a < k < kc, where k'^ = (1 -I- 2Qs), 
with A +00 as fc — > a -|- 0. On the other hand, A is 
negative and therefore no instability takes place at fc < a. 

Thus, this simplified analysis shows clearly that the 
instability found in Sec. IIV Al is generic. This kind of 
instability is not a feature of the particular problem, it 
is appears for any nontrivial distribution of the pulsation 
potential ?Ao- 



APPENDIX B: STABILITY OF 
QUASI-EQUILIBRIUM IN THE LIMIT OF 
LARGE NEGATIVE Qs 

To study the stability of the quasi-equilibrium state at 
large \Qs\ we use the method of matched expansions. We 
introduce a fast coordinate f = z/e. As before, = 
a/IQsI, for the sake of brevity we also suppress tilde for 
^. The solution of the inner problem depends on ^, and 
is sought in the form 



i) - 



i)a + e(j)i + e^(j)2 + 
(^«+^«+...J, 



(Bl) 
(B2) 

(B3) 



The solution of the outer problem, which depends on z, 
is presented as 



(o) _ 



e.s.t 

e 



( (o) 



(B4) 
.(B5) 



Here, ^'e.s.t." is used to denote exponentially small terms. 
Since 0*-°-' is negligibly small, we omit the superscripts for 
, j=0,l,2,.... 
Next, we assume that the growth rate is large in the 
sense 



A = e-'A 



(B6) 



and also take into account power expansions of $o and 
■00 given by relations ((2T|) with respect to e. As a result 
we obtain 



- A0O 


- 0, 


(B7a) 




= 0, 


(B7b) 




= 0, 


(B7c) 
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where subscript ^ is applied to denote the derivative with 
respect to f . 

By means of an ansatz 0o = 0o exp , Eq. (jB7a| 

is reduced to Hermite's equation: 



(B8) 



which for A„ = —2n, n = 0, 1, 2, . . . admits the solution 
given by the Hermite polynomials. Other possible values 
of A and the corresponding solutions are out of interest 
because no proper matching with the outer problem can 
be achieved. 

Accounting for the rescaling of the growth rate, see 
relation (|B6|) . we end up with result (|^^ for the spectrum 
of growth rates. The solutions with n > describe very 
fast temporal decay of perturbations. Hence, the only 
case that should be analyzed separately corresponds to 
n Ao = 0, when A = 0(1). In this case, the solution of 
Eq. (jB7a[) is as follows 



(B9) 



so that 00 coincides with <i>o"'', cf. Eq (|21ap . 
Solutions of Eqs. (jB7bp and (|B7cp are given by 



(BIO) 



where Bi and B3 are constants and g{(,) is as in Eq. (|21cp . 
Note that because of symmetry the quadratic and con- 
stant terms with respect to ^ are vanishing in the solution 
for the streamfunction. 

To the first order we obtain 

0145 + 2(^01)^ = -2[4°^(^*o)5 + (^5)5'/'o)B,lla) 



j4 (i) 



6ifc$aC0o, 



(Bllb) 



and to the second order we arrive at 



02« + 2 (C02)5 = + (A + fc2) 00 + zA;<^o$[,JC,(B12) 

where F is the term unimportant for the further analysis. 
This term includes the first order corrections to 4'*^*' and 
the second order corrections to the base state. The first 
order correction to the potential as well as the second 
order of the streamfunction are not needed below. 

The solvability condition for Eq. (jBlla[) can be ob- 
tained by integration of the equation over ^ from zero 
to infinity. Thus, Eq. (|Blla|) is solvable. However, its 
solution is not used below and for this reason is not pro- 
vided here. A similar condition for Eq. (|B12p leads to a 
relation 



A -I- A:^ - ikBi + ikBg 



0. 



(B13) 



The constants Bi and B3 entering Eq. (|B13|) should be 
found by means of the matching procedure. The correc- 
tion to the streamfunction is given by 



^d77 r giOdC 
Ja 



(B14) 



Keeping in mind the behavior of g(^) at large ^, one 
obtains an asymptotical law 



(0 



(B15) 



Hence, the solution of the inner problem for the stream- 
function at large ^ is: 



(B16) 



This solution must be matched with the solution of the 
outer problem: 



(B17) 



with the no-slip condition at z = 1. Since the perturba- 
tions of the concentration are exponentially small in the 
bulk, no external force acts on the fluid in this domain. 
The solution of Eq. (|B17p that satisfies the boundary 
conditions at z = 1 is 



Ci 



/ sinh/czi sinhgzi 



C2 (cosh kzi 



q 

coshqzi) 



(B18) 



where 21 = 1 — z and = k"^ + XS^^. By expanding 
this solution near z = and equating the coefficients at 
equal powers of z with those in Eq. (IB16p . we find that 



Bi = 3ifc$a 



with Qj. = 



sinh k sinh q — 2kq (cosh k cosh g — 1) 
q^ {q cosh q sinh k — k cosh k sinh q) 



(B19) 



0. 



Bearing in mind that q and q± depend on A, we sub- 
stitute these constants into Eq. (jBlSp and obtain a tran- 
scendent equation with respect to A. In the practically 
relevant case of S* » 1, it is necessary to expand q near 
k, which results in Eq. Note that this approxima- 

tion works well already at = 100, see Fig. [SJ^b). More 
precisely, the line corresponding to formula (j45p cannot 
be distinguished from the numerical results based on the 
solution of Eqs. (|Bl3l) and (|Bl9l) . 
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